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Existence, stability and dynamics of soliton complexes, centered at the site of a single transverse 

link connecting two parallel 2D (two-dimensional) lattices, are investigated. The system with the 

on-site cubic self-focusing nonlinearity is modeled by the pair of discrete nonlinear Schrodinger 

equations linearly coupled at the single site. Symmetric, antisymmetric and asymmetric complexes 

are constructed by means of the variational approximation (VA) and numerical methods. The 

VA demonstrates that the antisymmetric soliton complexes exist in the entire parameter space, 

while the symmetric and asymmetric modes can be found below a critical value of the coupling 

parameter. Numerical results confirm these predictions. The symmetric complexes are destabilized 

(^ ' via a supercritical symmetry-breaking pitchfork bifurcation, which gives rise to stable asymmetric 

CNj , modes. The antisymmetric complexes are subject to oscillatory and exponentially instabilities in 

narrow parametric regions. In bistability areas, stable antisymmetric solitons coexist with either 

^ ' symmetric or asymmetric ones. 

>—i ; 

fT^ , PACS numbers: 03.75.Lm; 05.45.Yv 



00 



> 



I. INTRODUCTION 



Solitons trapped at interfaces between different nonlinear media [l|, [4I, or pinned by defects [S-Q; have been the 
subject of many recent studies. It has been found that the self-trapped surface modes possess novel properties in 
comparison with the solitons in bulk media. Among noteworthy features of these localized modes are a threshold 
value of the norm, above which they exist, and the coexistence of different surface modes with equal norms. 

The studies of interface solitons in discrete systems have shown that spatially localized states with broken symmetries 
can exist [5|, which may be related to the general phenomenon of the spontaneous symmetry breaking (SSB) in bimodal 
nonlinear symmetric settings with a linear coupling between two subsystems. For the first time, the SSB bifurcation, 
which destabilizes symmetric states and gives rise to asymmetric ones, was predicted in a discrete self-trapping model 
^/~^ , in Ref . [9|] . This finding was followed by the prediction of the SSB in the model of dual-core nonlinear optical fibers 
CsJ ' [l3l) UM- ^^ the framework of the nonlinear Schrodinger (NLS) equation, the concept of the SSB was, as a matter of 

1^ [ ' fact, first put forward in early work 12]. Related to this context is the analysis of the SSB of discrete solitons in the 

f— V ■ system of linearly coupled ID and 2D (one- and two-dimensional) discrete nonlinear-Schrodinger (DNLS) equations 



[13j (the general outline of the topic of DNLS equations was given in book |14| ' 

T— I , The SSB was also analyzed in detail for solitons in the continual model of dual-core fibers with the cubic (Kerr) 

ILJ ' nonlinearity }15]-11.7|i and in related models of Bose-Einstein condensates (BEC's) loaded into a pair of parallel-coupled 

.^H [ cigar-shaped traps [T^l. In the latter context, the analysis was generalized for 2D coupled systems [I^. The SSB for 

J^ ■ gap solitons was studied too, in the model of dual-core fiber Bragg gratings [20|, and later in the model of the BEG 

J-j \ trapped in the dual-trough potential structure, combined with a longitudinal periodic potential [21|. As concerns the 

. . . 1 relation between discrete and continual systems, it is relevant to mention work [22l |. where exact analytical solutions 

were found for the SSB in the model with the nonlinear coefficient in the form of a pair of delta-functions embedded 

into a linear medium. In its own turn, the latter model has its own discrete counterparts, in the form of a pair of 

nonlinear sites embedded into a linear chain [2J, [2J| , or side-coupled to its \2U\ . These models may be realized in 

terms of BEG and optics alike. The SSB in such settings was recently analyzed in Refs. [2J| and ;25||, respectively. 

Goming back to discrete media, the objective of the present work is to study localized modes at the interface of 
two 2D uniform lattices with the cubic on-site self-focusing, which are linearly linked at a single site. This link plays 
the role of the interface. Accordingly, the localized modes are complexes formed by two fundamental solitons in 
each lattice, centered at the linkage site. We focus on the symmetry of the soliton complexes, with the intention to 
investigate the SSB transitions in them. 

This work is a natural extension of the recent study of interface modes in the system of single-site-couplcd ID 
nonlinear lattices [26]. Unlike the ID situation, it would be difficult to realize such a 2D setting in optics. However, it 
is quite possible in BEG: one may consider two parallel pancake-shaped traps combined with a deep 2D optical lattice 
traversing both pancakes [23, with the local link induced by a perpendicular narrow laser beam. In fact, similar 
two-tier layers of nonlinear oscillators, transversely linked at sparse sites, can be realized in a number of artificially 
built systems. 



The article is organized as follows. The model is formulated in Section 2. In the same section the existence and 
stability of various on-site-centered fundamental localized modes are considered in a quasi-analytical form by means of 
the variational approximation (VA) and Vakhitov-Kolokolov (VK) stability criterion. Numerical results for the soliton 
complexes are presented in Section 3, including the stability and dynamics. The numerical findings are compared to 
the predictions of the VA, and both the analytical and numerical results are compared to those reported in Ref. [2^] 
for the fundamental localized modes in the ID counterpart of the system with parallel-coupled lattices, or 'system 2', 
in terms of Ref. [2a]). The paper is concluded by Section 4. 



II. THE MODEL AND VARIATIONAL APPROXIMATION 



The formulation 



The set of the locally linked 2D uniform lattices is displayed in Fig. [T] The intra-site coupling constant in the 
lattices is C > 0, and e > is the strength of the transverse link. The lattice system is modeled by the following 
DNLS system. 



C , 



dt 



L'n+ 



l,m + (f>n-l,m + (f>n,m+l + 't'n.m-l 



+ e'0n,m<^n,O^m,O +7l0«,m| 4>n 



0, 



i 4^ + -^(V'n-(-l,m + V'n-l,m + ^n,m+l +^n,m-l - 4VAi,m) +£0„,m<5„,O'5m,O +7|'0„,^ '0n,m = 0, (1) 

where 7 > is the coefficient of the on-site self-focusing nonlinearity, t is the time, and Sm,n is the Kronecker's symbol. 
By means of obvious rescaling, we set C/2 = 7=1. 
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FIG. 1: (Color online) A schematic presentation of parallel 2D identical latices, linearly linked at the single site n — m = 0, 
with the coupling constant e. 



To construct soliton complexes formed around the transverse link, we look for stationary solutions, 4'n,m = 
Un.mGxp{—ifit), Tpn.m — Vn,m sxp {—i^t) , wherc Un,m,Vn,m and jjL are real lattice fields and the propagation con- 
stant, respectively. The corresponding stationary equations following from Eqs. ([1]) are 






Equations ^ can be derived from the Lagrangian, 
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where L„ and it, are the intrinsic Lagrangians of the uncoupled lattices, and the last term in Eq. ([3]) accounts for 
the coupling between them. 

B. The variational approximation 

The variational method follows the route described in Ref. [ij|. We adopt a natural ansatz, which was first applied 
to discrete lattices in Ref. [2^: 

{um,n,Vn^rn} = {A, B} cxp (-a|n|) cxp {-a\m\), (5) 

where A and B (but not a, see below) are treated as variational parameters. This form of the trial function admits 
different amplitudes, A ^ B, oi the solutions in the coupled lattices, but postulates equal widths in both of them, 
a~^. In this context, the SSB is signaled by the emergence of asymmetric solutions, with A^ ^ B^ p^[20|. 

The inverse width a of the localized trial solution is found independently from the linearization of Eqs. ([2]), which 
is valid in the soliton's tails (at |m|, \n\ — > oo), 

a = - In ((4 - m)/4 - ^(4 - lif /IQ - l) , (6) 

provided that the propagation constant is negative, /i < (otherwise, the solution cannot be localized). Relation (j6]) 
may also be cast in another form, that will be used below: 



.-e-"^i^-^(i^-l,, = 4-2(. + .-). (7) 

The substitution of ansatz ^ into Eqs. ([3]) and (|H) yields the corresponding effective Lagrangian with two varia- 
tional parameters A and B, where Eq. ([7]) is used to eliminate /i in favor of s: 

icff = (i«)eff + (i«)cff + 2£^S, (8) 






2 (1 - .s4)2 



{L,U = -2B'— + -B^ ' '' (9) 



The Euler-Lagrange equations for amplitudes A and B are [d/dA) [L^)^^ + 2eB = 0, {d/dB) {Li,)^^^ + 2eA = 0, or, 
in the explicit form, 



l + s\ , (l + s^) 



3 



-2 yl+^- j^^rA-'+eB = 



2- 



S (1-S4) 

-2^^ + 11^^^ + -^ = 0- (10) 

These equations allow us to predict the existence of three different types of the complexes formed by the fundamen- 
tal localized modes centered at the linkage site: symmetric and antisymmetric ones, with A = B and A = —B, 
respectively, and asymmetric modes with A^ ^ B^. 

1. Existence regions for the interface soliton complexes 
The solution for the symmetric soliton complexes is easily obtained from Eqs. ([T0| : 



^.^d-' 



4\2 



(1 



-(l + s2)-e 



(11) 



PQ 

< 




FIG. 2: (Color online) Amplitudes A and B of asymmetric, symmetric, and antisymmetric solitons (red, black, and blue colors, 
respectively), as predicted by the variational approximation (lines) and obtained in the numerical form (triangles, circles, and 
squares pertain to the asymmetric, symmetric, and antisymmetric modes, respectively) vs. the inter-lattice linkage strength e. 
The propagation constant is fixed to /x = —5. The dotted green vertical lines denote the numerically found critical values of 
e limiting the existence regions of the asymmetric (ec) and symmetric (se) solitons. Filled and empty symbols correspond to 
unstable and stable solitons, respectively, as concluded from the numerical investigation. 



As follows from Eq. ([TT|) . the existence domain of the symmetric solutions is e <£e = 2(1 + s'^)/s. For /i = —5, 
solution branches of all the types, produced by the VA along with their numerical counterparts, are shown in Fig. [51 
and the respective existence regions, including the one given by curve Se{p), is displayed in Fig. [3l The procedure for 
obtaining numerical results is described below. 

The amplitudes and existence range of the asymmetric solution complexes, with A^ ^ B^, can be calculated by 
adding and subtracting Eqs. (jlOp . After a straightforward procedure, the following expressions for the amplitudes 
are obtained: 



A 
B 



= ± 



rT^v/i + -^^ + v/aT^ 



(1 - s^ 



(1 



c4^2 



A 



(12) 



This soliton mode exists at £ < ^c = (1 + s^)/s, where Ec is the bifurcation value. The existence of the asymmetric 
mode may be naturally expected when the linkage constant (e) is not too large [ij], [1^1- [3 • Indeed, the ultimate 
asymmetric solution, with A ^ and i? = 0, is obviously possible in the limit of e = 0, which corresponds for the 
decoupled lattices. It is also natural that, with the decrease of e, the symmetry-breaking bifurcation should occur at 
some e = Sc, where the two asymmetric branches emerge from the symmetric one, see Fig. [21 However, the stability 
of the related solutions cannot be predicted solely by the VA. 



15 



lo- 



co 







1 








. 


□,' 






\ 


- 


SyS stable 


□'- 


8 

•- e 


•_ 




" AS does not exist 


~^-, 




~ 


■--o_ 




E 


°'"V' 






~^>- 


c 




•, 


SyS unstable 






~^>-^ 


° "''-, 


AS stable 








^^'^^'^^--^ 


' 1 








■ ' 



-8 







1^ 



FIG. 3: (Color online) The existence and stability diagrams for the fundamental symmetric (SyS), asymmetric (AS) and 
antisymmetric (AnS) solitons. The variational results are presented by curves, and numerical findings — by symbols in the 
parameter space {e,n). Black circles and the corresponding line mark the boundary of the AS existence region. Red squares 
and the dashed line denote the existence boundary for the SyS modes. The numerical calculations show that the symmetry- 
breaking bifurcation takes place along the curve £c(m). The AnS mode exists in the entire parameter plane. According to the 
Vakhitov-Kolokolov criterion, the AS may be stable in the whole existence region, the SyS may be stable above the orange 
dashed line, and AnS — to the left from the blue dotted line. The stability type of the SyS and AS modes, indicated in the 
figure, was established by means of numerical computations. 



For antisymmetric soliton complexes, the amplitude is obtained from Eq. (|10p by setting A 



A' = 



{i + s^y 



s 



-B, which yields 

(13) 



cf. Eq. (jlip for the symmetric modes. Relation ([T3| is plotted versus e for fixed /j, = —5 by the blue (upper) curve 
in Fig. [21 As follows from this relation, the VA predicts the existence of the antisymmetric solitons in the entire 
parameter space, on the contrary to the limited existence regions predicted for the symmetric and asymmetric modes. 



2. Stability of the soliton complexes 



The stability of the discrete solitons predicted by the VA can be estimated by dint of the VK criterion, dP/dfi > 

^ ^) is the total norm (power) of the soliton complex [23|. The norm 



Q, ^\^e.e P ^ YZT-ZTZT-Ziul.rn 
corresponding to the ansatz ([5]) is 



P = {A^ 



B' 



For the symmetric solution with A = B, Eqs. (|11|) and ([T4| yield 



P = 2 



(1 + 5^ 
(1 + S4 



s 



(14) 



(15) 



This expression satisfies condition dP/ds < 0, which is tantamount to the VK criterion, in the region of 

.>^(i±i!l-ii±^!lii±i!l. (16) 

This region is displayed in Fig. ^ as the area above the dashed orange curve (in the right bottom corner of the 
figure). Thus, according to the VK criterion, the stable symmetric branch may exist in the large part of the parameter 
plane, except for the small region below the dashed orange curve, which corresponds to the weak inter-lattice linkage 
and wide solitons. 

For the asymmetric solutions, the substitution of Eq. p^ into Eq. p^ produces a simple expression for the total 
norm, which does not depend on e: 

(1 + s2)5 

P-^ \,\ L - (17) 

s(l + S^)"^ 

It also satisfies the VK criterion in the entire region of the existence of the asymmetric mode. 
For the antisymmetric modes with A ~ —B, the use of Eq. ((T5| gives 



P. 2-^1 



(1 
In this case, condition dP/ds < is satisfied at 



2, 9x 



(18) 



.<(i±i!l^±i!)-^ii±i!l. (19) 

■is-^ s 

The region defined by Eq. p9|) is shown in Fig. ([3]) as the area to the left from the blue dotted curve (the one which 
cuts the entire plane). 

However, the VK criterion offers only the necessary condition for the soliton stability. To identify regions of full 
stability of the soliton complexes, the VK criterion should be combined with the spectral condition, which requires 
the existence of only pure imaginary eigenvalues in the linearization of Eqs. ([T]) with respect to small perturbations 
around the stationary soliton solutions. The spectral analysis is performed numerically in the next section. In fact, 
the results demonstrate that, unlike the prediction of the shape of the soliton modes by means of the VA, the stability 
prediction based on the formal application of the VA criterion is not accurate. 

III. NUMERICAL RESULTS 

The predictions of the VA were verified by numerically solving stationary equations ([2]) . The numerical algorithm 
is based on the modified Powell minimization method [301 ■ The initial guess to construct soliton complexes centered 
at the linkage site of the 2D lattices (Fig. [Ij was taken a.s uq = vq = A > for symmetric solutions, uq = A > 
0, vq = B ^ A > for asymmetric ones, and uq = A > 0, vq = —A for solutions of the antisymmetric type, with the 
VA-predicted values of A and B. At other sites, the amplitudes of the initial guess are set to be zero. 

A. Stationary soliton modes 

The stability of the stationary modes was investigated through the calculation of eigenvalues (EVs) of small pertur- 
bations around the stationary solutions, which were computed following the lines of Refs. [2^, |30, |3l| . The obtained 
results were further verified in direct numerical simulations of Eqs. ([T]), using the sixth-order Runge-Kutta algorithm, 
as in Refs. [26 , 30]. The simulations were carried out for stationary soliton complexes to which initial perturbations 
were added. 

Typical shapes of symmetric, asymmetric and antisymmetric soliton complexes found in the numerical form are 
displayed in Fig. H) The respective dependencies of the amplitudes A and B on e were displayed above, for all 
the soliton branches, alongside their VA-predicted counterparts, in Figs. [2] and [3] The numerical results show that 
the symmetric and asymmetric complexes exist in bounded regions of the parameter space, see Fig. [31 The results 
predicted by means of the VA are in good agreement with their numerical counterparts for all the types of the soliton 
complexes. The comparison of the variational and numerical dependencies of the solitons' amplitudes on e, for fixed 
s ~ 0.23, which corresponds to /i = —5, can be seen in Fig. [2] The difference between the analytical and numerical 









FIG. 4: Typical shapes of the sohton complexes of different types: symmetric in (a) and (b), asymmetric in (c) and (d), 
and antisymmetric in (e) and (f). The parameters are /i = —3, e — 2.5, the corresponding norms for complexes being 



results is negligible for small e, and slightly grows with e. In particular, the existence border of the symmetric 
complexes predicted by the VA is eg = 9, while its numerical counterpart is Se = 8.63. 

Finally, it is easy to check that both the variational and numerical results demonstrate that the SSB pitchfork 
bifurcation observed in Fig. [2]is of the supercritical type, similarly to that reported in Ref. [26[ for the ID counterpart 
of the present system (in contrast with the subcritical bifurcation demonstrated by the system of ID and 2D parallel 
chains with the uniform linear coupling acting at each site ilSj). 



B. The stability analysis 



The linear-stability analysis demonstrates that the symmetric complexes emerge as stable solutions at e — Se, and, 
as expected, they change the stability at the bifurcation point, s — Sc, where the asymmetric solution branches appear 
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(see Figs. [2] and [3]). Unstable symmetric solutions are characterized by the pure real EV pairs, see Fig. [Sja). For the 

Mameter set used in Fig. [SJa), the isolated discrete solitons existing in the uncoupled lattices at e = are stable 
l32 | . The introduction of the linkage between the lattices changes the stability of the symmetric complex formed 
by such solitons. The results do not confirm the prediction, based on the VK criterion, that the symmetric complexes 
change the stability at the value of the e given by Eq. (TlG)) . which corresponds to the orange dashed curve in Fig. 
[3l The instantaneous destabilization of the symmetric bound states with the increase of e from zero [see Fig. El^a)] 
is simply explained by the fact that s ymm etric complexes are always unstable against the SSB in dual-core systems 
with a small linear-coupling constant 
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FIG. 5: (Color online) (a) Pure real eigenvalues (EVs ) vs. e for symmetric solitons with /i = —5. Dotted lines denote the 
bifurcation value Ec and the existence threshold e^ for the symmetric solitons. Pure real EVs, and the real part of complex 
EVs, are shown vs. e in panels (b) and (c), respectively, for antisymmetric solitons. Black (empty circles), orange (squares), 
and red (filled circles) symbols correspond, respectively, to fixed ^ — —1.5, —2, and —3. Dotted lines marked by £1,2,3 in plot 
(b), and by e' in plot (c) are boundaries of regions where the pure real EVs (b), or real parts of the complex EVs (c), take 
significant values, Re(EV) > 0.001. 

Direct simulations confirm the stability of the SyS complexes in the interval oi Sc < s < Se- On the other hand, 
simulations of the evolution of unstable symmetric modes demonstrate that, under the action of small perturbations, 
these unstable modes (at e < Ec) evolve into asymmetric breathing complexes, which consist of two oscillating localized 
components, that exchange energy in the course of the evolution. This behavior is illustrated in Fig. [51 where the 
evolution of the unstable SyS mode into the AS complex is shown by plotting the amplitudes of the component solitons 
versus time. 




FIG. 6: The evolution of the amplitudes of the unstable SyS mode with fi — —5, e — 3.4 into the corresponding stable AS 
complex. The amplitudes of the components of the latter complex are represented by different lines. 



Two mutually symmetric branches of asymmetric solutions, which are created by the destabilization of the sym- 
metric branch, see Eq. (|12p . turn out to be stable, according to the linear-stability analysis, which yields for them 
EV spectra with the zero real part. Direct simulations corroborate that slightly perturbed asymmetric complexes are 
robust modes, see Fig. [T] 

The purely real EV pairs which are numerically calculated for the antisymmetric modes become significant above 










FIG. 7: The oscillatory evolution of a perturbed asymmetric complex with P — 12.02, s — 5, and jj, = —8. The relative strength 
of small perturbation with respect to the solution amplitude is 0.01. Profiles of the two components are shown, in the top and 
bottom rows. 







FIG. 8: The time snapshots illustrating the evolution of the AnS complex with /x — —1.5, e — 7.2, and P — 33.48. Four plots 
represent the amplitude profiles of the component solitons, 10m, n| = IV"™."!) at i = 0,50, 100, and 200. The components are 
strongly pinned, with the amplitudes at the interface sites slightly oscillating around value 3.8. 



some threshold values of e, which depend on /i. The threshold values, £1,2,3 for /i = —1.5, —2, —3, are shown in Fig. 
^h) . On the other hand, the antisymmetric modes give rise to the complex EVs at £ < e' and arbitrary /j,, as shown 
in Fig. Oc). Therefore, at large values of the coupling constant, e > £1.2,3, the instabiUty is determined by the purely 
real EVs, but in the region of £ < min{e', £1.2,3) the real part of the complex EV dominates the instabiUty. Actually, 
the numerical results for the (in)stability of the antisymmetric complexes do not corroborate the analytical prediction 
presented by Eq. P^ . 

The dynamics of the antisymmetric complexes with positive real parts of the EVs is not significantly affected by small 
perturbations, except in the area with very small £, where the AnS, SyS and AS with close values of the power coexist, 
see Fig. [3] for e ^ Ec- In the latter case, the AnS complexes follow the same scenario as the unstable SyS complexes, 
i.e., the unstable AnS evolves into a breathing AS complex. The reason for the robustness of other antisymmetric 
modes is the fact that the corresponding branch features large amplitudes of the solitons, leading to their strong 
trapping at the central lattice sites. Therefore, the actual instability (escape of the discrete wave fields) is suppressed 
by the strong the Peierls-Nabarro potential barrier [33| . Accordingly, due to the weak effective coupling between the 
central site and the adjacent ones, the introduction of small perturbations can excite internal oscillations but does 
not destroy the localization of the mode, and the emerging breather (which keeps the antisymmetric structure, as 
concerns the relation between its components) remains a strongly trapped state. This is the case in almost the whole 
existence region of the antisymmetric solitons. The dynamics of exponentially unstable antisymmetric complexes 
(actually, with large pure real EVs) is illustrated in Fig. [S]for /i = —1.5 and e = 7.2, with norm P = 33.48. Keeping 
the antisymmetric structure, as said above, the two components feature identical amplitude profiles, |V'n,m| — \4'n,m\- 
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Both of them shrink into more pinned modes that radiate away a small part of their norm (energy) , which forms a 
small but finite oscillating background, while the central peak remains robust. 

The instability of the antisymmetric modes in the case of very small pure real EVs develops very slowly, in the 
presence of small perturbations, following the same scenario. In Fig. 1^1 the evolution of a typical antisymmetric 
complex subject to the oscillatory unstable is shown. The part of the initial energy lost into background is smaller 
then in the previous case, and the central peak slightly oscillates in time. 
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FIG. 9: The evolution of a perturbed antisymmetric mode which is unstable against oscillatory perturbations, with fi — —3, 
e = 1.5, and P — 17.64. The component solitons are strongly pinned, and their amplitudes slightly oscillate, similar to the case 
displayed in Fig. [8] 



Returning to the global existence diagrams, it is worth to note that two bistability areas can be identified in 
them: the domain of the coexistence of stable symmetric and quasi-stable antisymmetric solitons (the quasi-stability 
pertains to very small growth rates mentioned above), or the one featuring the simultaneous stability of asymmetric 
and antisymmetric modes (at very small e), on the opposite side on the SSB bifurcation. This result is in accordance 
with similar findings reported in other linearly-coupled two-component systems with the self-focusing nonlinearity 

Eliilii. 

It is also relevant to compare properties of the soliton complexes in the present two-component 2D lattice system, 
and on-site solitons in the uniform 2D lattice described by the single DNLS equation, which corresponds to Eqs. ([2]) 
with e = [30, ll^l- Along the line of e = in Figs. [2l [3] and[5l which correspond to fi — —5, one finds a stable 
symmetric complex, a stable asymmetric mode (with one zero component), and stable antisymmetric complexes. In 
terms of the uniform 2D lattice, the symmetric complex is formed of two identical on-site-centered discrete solitons, 
which are stable in the usual DNLS lattice [30, ^^ [see, e.g.. Fig. 5(a) in Ref. [30]]. As mentioned above, the 
introduction of the intcr-latticc linkage (e > 0) leads to the onset of the exponential instability in the complex formed 
by two identical fundamental on-site solitons. The symmetric complex recovers its stability at s > Ec- Therefore, one 
can associate two bifurcation points with values e = and e = Sc- The latter one was actually identified above as 
the supercritical pitchfork bifurcation at which the stable asymmetric branches disappear and the symmetric one is 
re-stabilized. 



IV. CONCLUSION 



In this work, we have introduced the 2D double-lattice nonlinear system, linked in the transverse direction at a single 
site. The system can be realized in terms of BEC, and may actually occur in a range of artificially built discrete non- 
linear media. The on-site nonlinearity was considered to be self-focusing (the self-defocusing can be easily transformed 
the same system by means of the staggering transformation IJ]). We have used the VA (variational approximation) 
and numerical methods to find the regions of existence, in the parameter plane of (/x, e) (the propagation constant 
and the strength of the linkage between the lattices), of the localized symmetric, asymmetric and antisymmetric 
soliton complexes pinned to the linkage site. It was shown, by means of both approaches, that the existence regions 
of the symmetric and asymmetric complexes are bounded. The SSB (spontaneously symmetry-breaking) pitchfork 
bifurcation of the supercritical part has been found, which destabilizes the symmetric complexes and simultaneously 
creates stable asymmetric ones. The stability of the antisymmetric complexes changes twice. Areas of the bistability 
between the antisymmetric modes and either symmetric or asymmetric ones have been found too. Direct simulations 
demonstrate that unstable symmetric modes relax into the breathing asymmetric complexes, while the antisymmetric 
solitons with large amplitudes are robust against perturbations, transforming into strongly pinned breathers, which 
keep the antisymmetric structure. 
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This work also suggests a possibility to create and investigate vortex complexes in 2D parallel-coupled lattice 
systems, "which will be reported elsewhere. 
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